%% The Discrete Fourier Transform
%
%
% Any vector $X(n)$ can be represented by the finite sum:

%%
%
% $$ X(n)=\frac{1}{N}\sum_{k=1}^N Y(k) e^{2\pi i (n-1)(k-1)/N}$$
%
% where the coefficients $Y(k)$ are given by the Discrete Fourier Transform (DFT):
%%
%
% $$Y(k)=\sum_{n=1}^N X(n) e^{-2\pi i (n-1)(k-1)/N}$$
%
% Note the $k$ used here means the 'spatial frequency', or
% $k=\frac{1}{\lambda}$
% and therefore the $2\pi$ in the exponent is required.
% Spacing in frequency domain is given by 1/L, where L is the length of the vector.


%% Calculating the Coefficients
%
% We can see that each coefficient Y(k) is the sum of the element-wise
% product of the vector $X(n)$ and a the vector $e^{-2\pi i (n-1)(k-1)/N}$
%
% An animation of the calculation of $Y(k)$

Lambda=10;
N=50;
SamplingRate=1;
NVec=(1:N)*SamplingRate;
X = sin(2*pi/Lambda*(NVec-1));
h=figure
Y = fft(X);

for kk=1:N
    M=exp(-2*pi*i*(NVec-1)*(kk-1)/N);
    
    subplot(2,2,1)
    hold off
    plot(NVec,X,'b')
    hold on
    plot(NVec,real(M),'r')
    xlabel('n')
    ylabel('X(n)')
    legend('X','real(M)')
    
    subplot(2,2,2)
    hold off
    plot(NVec,X,'b')
    hold on
    plot(NVec,imag(M))
    xlabel('n')
    ylabel('X(n)')
    legend('X','imag(M)')
    
    subplot(2,2,3)
    plot(1:kk,real(Y(1:kk)),'r')
    xlabel('k')
    ylabel('Y(k)')
    legend('real(Y)')
    axis([1 N min(real(Y)) max(real(Y))])
    
    subplot(2,2,4)
    plot(1:kk,imag(Y(1:kk)),'g')
    xlabel('k')
    ylabel('Y(k)')
    legend('imag(Y)')
    axis([1 N min(imag(Y)) max(imag(Y))])
     
    %pause(1)
end

%% Amplitude and Phase
%
% Remember $<sin(\omega t)^2>_t=<cos(\omega t)^2>_t=1/2$.  

TotalTime=100;
TimeStep=1;
N=TotalTime/TimeStep;
Tvec=(0:N-1)*TimeStep;

T=20; %Wave Period
F=1/T; %Wave frequency

Phase=pi/2
X= cos(2*pi*F*Tvec+Phase);
figure;plot(X)
xlabel('t')
ylabel('X(t)')

Y=fft(X)
K=(0:N-1)/TotalTime;
figure
plot(K,real(Y))
hold on
plot(K,imag(Y))
legend('real(Y)','imag(Y)')
xlabel('f')
ylabel('Y(f)')

Mag = abs(Y(6))
Phase = angle(Y(6))

%% Nyquist Frequency and Negative Frequencies
%
% Frequencies above half the sampling rate (Nyquist frequency) 
% can be seen as negative frequencies.  
%

TotalTime=10;
TimeStep=1;
N=TotalTime/TimeStep;
SamplingRate=1/TimeStep
Nyquist = SamplingRate/2

K=SamplingRate/N*(0:N-1)
K(K>Nyquist)=-1+K(K>Nyquist)


% An example showing the comparison

TotalTime=100;
TimeStep=1;
N=TotalTime/TimeStep;
Tvec=(1:N)*TimeStep;
SampleRate=1/TimeStep
Nyquist=SampleRate/2;

T=TotalTime/5; %Period of wave
F=1/T;
X = sin(2*pi*F*Tvec);
F2=-F;
X2 = sin(2*pi*F2*Tvec);
F3=SampleRate-F
X3 = sin(2*pi*F3*Tvec);

figure
plot(Tvec,X,'b')
hold on
plot(Tvec,X2,'go')
plot(Tvec,X3,'r','linewidth',2)
xlabel('t')
ylabel('X(t)')
legend('F','-F','SampleRate-F')

% The third case looks the same as the first two because of 'aliasing'
% web('https://en.wikipedia.org/wiki/Aliasing')
% If we plot it at higher sampling rate, we can see how this happens

TotalTime=100;
TimeStep=.1;
N=TotalTime/TimeStep;
Tvec2=(1:N)*TimeStep;

X3 = sin(2*pi*F3*Tvec2);

figure
hold on
plot(Tvec,X2,'go')
plot(Tvec2,X3,'r','linewidth',2)
xlabel('t')
ylabel('X(t)')
legend('-F','SampleRate-F')


%% Shifting FFT so Zero is at Center

close all

%Make a square pulse signal
N=100;
T=(0:N-1);
S=zeros(N,1);
S(T>30&T<50)=1;

figure
plot(T,S)
xlabel('T')
ylabel('Signal')
title('Signal')

%The Fourier Transoform
F=1/N*(0:N-1);
FTS=fft(S);
figure; hold on
plot(F,real(FTS))
plot(F,imag(FTS))
xlabel('F')
ylabel('FFT')
legend('real','imag')
title('Un-Shifted Transform')

% Shift zero to center
FTS_Shift=fftshift(FTS);

F_Shift=1/N*(-N/2:N/2-1)% For even number of points
%F_Shift=1/N*(-floor(N/2):floor(N/2))% For odd number of points

%The shifted Transform
figure
plot(F_Shift,real(FTS_Shift))
hold on
plot(F_Shift,imag(FTS_Shift))
xlabel('F')
ylabel('FFT')
legend('real','imag')
title('Shifted Transform')


%% Fourier-Domain Filtering   
% Let's remove high frequencies

Ncut=20;
FTS_Shift(1:Ncut)=0;
FTS_Shift(end-Ncut+1:end)=0;

% The Filtered, shifted Transform
figure
plot(F_Shift,real(FTS_Shift))
hold on
plot(F_Shift,imag(FTS_Shift))
xlabel('F')
ylabel('FFT')
legend('real','imag')
title('Filtered, Shifted Transform')

% Inverse transform
FTS_Filter=ifftshift(FTS_Shift);
S_Filter=ifft(FTS_Filter);

% Plot our filtered signal
figure
plot(T,S_Filter)
xlabel('T')
ylabel('Signal')
title('Filtered Signal')

%notice Gibbs phenomenon
web('http://en.wikipedia.org/wiki/Gibbs_phenomenon')










